Magnetic Anomaly Maps

This file is best viewed in a Pluto notebook. To run it this way, from the MagNav.jl directory, do:

julia> using Pluto
julia> Pluto.run() # select & open notebook

This is a reactive notebook, so feel free to change any parameters.

Import packages and DataFrames

The DataFrames listed below provide useful information about the flight data (collected by Sander Geophysics Ltd. (SGL) in 2020) & magnetic anomaly maps.

Dataframe	Description
`df_map`	map files relevant for SGL flights
`df_cal`	SGL calibration flight lines
`df_flight`	SGL flight files
`df_all`	all flight lines
`df_nav`	all navigation-capable flight lines
`df_event`	pilot-recorded in-flight events

begin
    cd(@__DIR__)
    # uncomment line below to use local MagNav.jl (downloaded folder)
    # using Pkg; Pkg.activate("../"); Pkg.instantiate()
    using MagNav
    using CSV, DataFrames
    using Plots: plot, plot!
    using Random: seed!
    using Statistics: mean, median, std
    seed!(33); # for reproducibility
    include("dataframes_setup.jl"); # setup DataFrames
end;

Load Perth map

This is the Perth map (at 800 m) as provided by Sander Geophysics Ltd.

begin
    map_gxf = MagNav.ottawa_area_maps_gxf()*"/Perth_Mag.gxf"
    p_mapS_800 = map_gxf2h5(map_gxf,800)
    p1 = plot_map(p_mapS_800)
end

Display Perth map in Google Earth

Display the Perth map in Google Earth by uncommenting below to generate a KMZ file, then open in Google Earth.

# map2kmz(p_mapS_800,"Perth")

Overlaid on the Perth map are most of the flight lines used to generate the map, which were collected during Flight 1004 (see readme) & Flight 1005 (see readme). The full list of SGL flights is in df_flight.

begin
    xyz_1004 = get_XYZ(:Flt1004,df_flight;silent=true) # load flight data
    xyz_1005 = get_XYZ(:Flt1005,df_flight;silent=true) # load flight data
    lines = [4019,4018,4017,4016,4015,4012,4011,4010,4009,4008,4007,4004,
             4003,4002,4001,421,419,417,415,413,411,409,408,407,405,403,401]
    for line in lines
        xyz = line in xyz_1004.line ? xyz_1004 : xyz_1005
        ind = get_ind(xyz,line,df_nav) # get Boolean indices
        plot_path!(p1,xyz.traj,ind;show_plot=false,path_color=:black)
    end
    p1
end

Load previously processed maps

Eastern Ontario & NAMAD maps have been processed & saved previously. The full list of maps is in df_map.

begin # e_mask contains Boolean indices with "real" map data (not filled-in)
    e_mapS_395  = get_map(:Eastern_395,df_map) # load map data
    e_mask      = trues(size(e_mapS_395.map)) # todo: update with v1.2.1
end;

begin
    xx_lim = extrema(e_mapS_395.xx) .+ (-0.01,0.01)
    yy_lim = extrema(e_mapS_395.yy) .+ (-0.01,0.01)
    n_mapS_395 = upward_fft(map_trim(get_map(MagNav.namad),
                            xx_lim=xx_lim,yy_lim=yy_lim),e_mapS_395.alt)
end;

Plot all Ottawa area maps

3 maps are overlayed with 3 different color schemes.

begin
    clims = (-500,500)
    dpi   = 50
    p2 = plot_map(n_mapS_395;clims=clims,dpi=dpi,legend=false)     # 395 m
    plot_map!(p2,e_mapS_395;clims=clims,dpi=dpi,map_color=:magma) # 395 m
    plot_map!(p2,p_mapS_800;clims=clims,dpi=dpi,map_color=:gray)  # 800 m
end

Plot Eastern Ontario altitude map CDF

Most of the map data was collected at altitudes between 215 & 395 m.

begin # the Eastern drape map contains an additional drape (altitude) map
    e_mapS_drp = get_map(:Eastern_drape,df_map) # load map data
    e_alts  = e_mapS_drp.alt[e_mask]
    alt_avg = round(Int,mean(e_alts))
    alt_med = round(Int,median(e_alts))
    alt_val = 200:500
    alt_cdf = [sum(e_alts .< a) for a in alt_val] / sum(e_mask)
    p3 = plot(alt_val,alt_cdf,xlab="altitude [m]",ylab="fraction [-]",
              title="altitude map CDF",lab=false,dpi=200)
end

Plot Eastern Ontario altitude map

Minimal areas have map data collected at 395 m or higher (colored spots), so a level map can be generated at 395 m using almost entirely upward continuation & minimal downward continuation.

begin
    p4 = plot_map(e_mapS_395;dpi=dpi,legend=false,map_color=:gray)
    e_mapS_395_ = deepcopy(e_mapS_395)
    e_mapS_395_.map[e_mapS_drp.alt .<= 395] .= 0
    plot_map!(p4,e_mapS_395_;dpi=dpi)
end

The trend of the map values agree, as expected.

Create a 3D map

All of the functionality in MagNav.jl can use 3D maps, which contain multiple stacked 2D maps with constant altitude spacing between each map level. There are 2 ways to create a 3D map, the first of which is to upward/downward continue to multiple altitudes using the upward_fft function. In this case, a filled Perth map at 800 m is stacked with the same map upward continued to 810 m.

begin
    p_mapS3D_1 = upward_fft(p_mapS_800, p_mapS_800.alt .+ [0,10]) # 3D map
    println(p_mapS3D_1.alt) # map altitude levels
    p9 = plot_map(p_mapS3D_1) # plot single map level (default is lowest)
end

The second way to create a 3D map is to combine multiple 2D maps using the map_combine function. In this case, the Eastern map (395 m) & Perth map (800 m) are stacked into a 3D map.

begin
    p_mapS3D_2 = map_combine([e_mapS_395,p_mapS_800]; # 3D map
                             n_levels = 2,
                             dx       = MagNav.get_step(p_mapS_800.xx),
                             dy       = MagNav.get_step(p_mapS_800.yy),
                             xx_lim   = extrema(p_mapS_800.xx),
                             yy_lim   = extrema(p_mapS_800.yy));
    println(p_mapS3D_2.alt) # map altitude levels
    p10 = plot_map(p_mapS3D_2) # plot single map level (default is lowest)
end